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1. Introduction 



To clarify the phase diagram of QCD and thus the nature of matter under extreme conditions 
is one of the most interesting and fundamental tasks of high energy physics. Lattice QCD has been 
shown to provide important and reliable information from first principals on QCD at zero density. 
However, Lattice QCD at finite densities has been harmed by the complex action problem ever 
since its inception. For jj. > the determinant of the fermion matrix (detM) becomes complex. 
Standard Monte Carlo techniques using importance sampling are thus no longer applicable when 
calculating observables in the grand canonical ensemble according to the partition function 

Zgc(m) = y"^f/detM[U](M)exp{-S G [U]}. (1.1) 

Recently many different methods have been developed to cirumvent the complex action problem 
for small fi/T [jl|, §]. For a recent overview see also [§]. 



2. Formulation of the method 



A very general formulation of the DOS method is the following: One exposed parameter (<j>) 
is fixed. The expectation value of a thermodynamic observable (O), according to the usual grand 
canonical partition function ( |l.l[ ), can be recovered by the integral 

<0>=jd<i> (0/(t/)> p(0) I Jd<$> </(I/)) p(0) (2.1) 

where the density of states (p) is given by the constrained partition function: 

p(x)=Z^x) = JmJg(U)S^-x). (2.2) 

With ( )^ we denote the expectation value with respect to the constrained partition function. In 
addition, the product of the weight functions f,g has to give the correct measure of Zee'- fg = 
detMexp{— Sq}- This idea of reordering the partition functions is rather old and was used in 
many different cases [Q, ||, ^| The advantages of this additional integration becomes clear, when 
choosing <p = P and g(U) = 1. In this case p(0) is independent of all simulation parameters. The 
observable can be calculated as a function of all values of the lattice coupling p\ If one has stored 
all eigenvalues of the fermion matrix for all configurations, the observable can also be calculated 
as a function of quark mass (m) and number of flavors^ (N/). In this work we chose 

0=P and g= |detM|exp{-5 G }, / = exp{j0}. (2.3) 

In other words we constrain the plaquette and perform simulations with measure g. In practice, 



we replace the delta function in Equation (2.2) by a sharply peaked potential HQ]. The constrained 



partition function for fixed values of the plaquette expectation value can then be written as 

p(x)aJmjg(U)txp{-V(?c)}, (2.4) 
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where exp{— V(x)} is a Gaussian potential with 

V(x)= l -y(x-Pf. (2.5) 

We obtain the density of states (p(x)) by the fluctuations of the actual plaquette P around the 
constraint value x. The fluctuation dissipation theorem gives 

^Anp(x)=<x-P> x . (2.6) 



Before performing the integrals in Equation ( |2T| ) we compute from an ensemble generated at 
(Mo,j3o): 

(Of(U)) x (n,f3) = (Of{U)R(ji,to,P,h))J(R(ji,i4>,p,Po)) x , (2.7) 
(f(U)) x (n,P) = (f(U)R(^^,P,Po)) x /(R(^^,P > p )) x , (2.8) 

^\np(x,n,P) = ((x-P)R(n,^,P,Po)) x . (2.9) 

Here R is given by the quotient of the measure g at the point (jtt , j3 ) and at the simulation point 

(jk>,Po), 

R(^^,P,Po)=g(H,P)/g(Ho,Po) = -M^l exp{5 G (j3)-5 G (j3 )}. (2.10) 

|det(jU )| 

Heaving calculated the expressions (p!7|)-( ^9| ), we are able to extrapolate the expectation value 



of the observable (2.1) to any point (jU,/3) in a small region around the simulation point (/io,j3o)- 
For any evaluation of (O) (£i,j3), we numerically perform the integrals in Equation ( |2.1[ ). We also 
combine the data from several simulation points to interpolate between them. 

3. Simulations with constrained plaquette 

The value we want to constrain is the expectation value of the global plaquette, which is given 
on every gauge configuration by the sum over all lattice points (y) and directions (/iv) of the local 
plaquette P^ v (y) and its adjoint Pji V (y), 

P = L I kiwVv(y)+TrfJv(y)]- 

y K(i<v<4° 

Since the plaquette is also the main part of the gauge action, 

5 c = -j3£ I U[Tr^ v W + Tr^ v W]-l), (3.2) 

x l<fi<v<4 I D J 

the additional potential V can be easily introduced in the hybrid Monte Carlo update procedure of 
the hybrid-R algorithm After calculating the equation of motion for the link variables Un(y), 
we find for the gauge part of the force 



P U,(y)T,(y)(l + ] 



(3.3) 

TA 



.3 MW MW V P 
Here the subscript TA indicates the traceless anti-Hermitian part of the matrix. We see that in 
each molecular dynamical step the measurement of the plaquette is required. However, the only 
modification in the gauge force is the factor in round brackets. 
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Figure 1: Results for Simulations at /3 = 4.98, fx = 0.3, A = 0.02, «/ = 4, am = 0.05, and number of lattice 
points: 4 4 . Shown are: (a) the density of states p(x), the phase factor (cos(0)}, and their product, (b) the 
Plaquette as a function of the coupling /3, with and without the phase factor, (c) the Susceptibility of the 
Plaquette as a function of the coupling j8, with and without the phase factor, and (d) the extrapolation J3 C (A) 
to A = 0, with and without the phase factor. 



4. The critical line and the determination of a triple-point 

Simulations have been performed with staggered fermions and Nf = 4. We chose 9 differed 
points in the (j3 , ju) -plane for the 4 4 lattice and 8 points for the 6 4 lattice. On each of these points we 
did simulations with 20-40 constrained plaquette values, all with quark mass am = 0.05. Further 
simulations has been done with (j3,ju) = (5.1,0.3) on the 6 3 x 8 lattice for am = 0.05 and am = 
0.03. 

Fist of all we check, whether we can reproduce old results with our new method. We show 
in Figure [l|(a) results from a Simulation at /I = 0.3, /3 = 4.98 and A = 0.02. Here we plot the 
density of states (p) and the real part of the phase factor {cos (6)) as a function of the constrained 
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Figure 2: The phase diagram in physical (a), and the quark number density at constant temperature T = 
143 MeV (4 4 lattice), T = 124 MeV (6 4 lattice) and T =93 MeV (6 3 x 8 lattice). 

plaquette value. The results have been interpolated in P, in order to obtain a better result for the 
necessary integration over P. The distribution p shows a clear double peak structure, which signals 
the transition. The phase factor is smaller in the low temperature phase (P<2.8). Hence in the 
product p (cos(d)) the low temperature peak is suppressed. Now we perform the integrals 

(P) = Jdxxp(x)(cos(e)) x , (P 2 ) = j dxx 2 p(x)(cos(d)) x . (4.1) 

In Figure [l|(b) we plot the plaquette expectation value (P) as a function of the coupling j3 . The 



j3-dependence is given by Equations (|27^-([L9|). We indeed find that including the phase factor 
does shift the transition to lower values of the coupling, which also means to lower temperatures. 
This can also be seen in a shift of the peak of the susceptibility of the plaquette %p = (f 2 ) ~~ (P) » 
which we plot in Figure [|(c). Since the A parameter introduces a systematic error, which can be 
seen by the relative large critical coupling of jS c = 4.976(4) in comparison to the result form multi- 
parameter reweighting j3 f = 4.938(2) [[]]], we perform the a linear extrapolation of A — ► 0, from 
A = 0.02, A = 0.015 and A = 0.01. We show the extrapolation in Figure |l|(d). The extrapolated 
result jS = 4.938(4) (including the phase factor) and the result from multi -parameter reweighting 
are in very good agreement. From now on we only give results for Xfm = 0.2, the A dependence 
is however expected to be smaller for larger }i. 

In the range of 0.4<aju<0.5 for the 4 4 lattice, as well as 0.3<a/i<0.4 for the 6 4 lattice we 
found two transitions in the plaquette expectation value (P) (j8). The two critical couplings result 
in two transition lines in the phase diagram. The two transition lines are almost perpendicular in the 
(/3,jU)-diagram, and join in a triple-point of the phase diagram, in Figure |2|(a) we show the phase 
diagram in physical units. The scale was set by the Sommer radius ro, measured on a 10 3 x 20 
lattice. In both cases, the triple point is located around w 300 MeV, however its temperature 
(r tri ) decreases from T tri ss 148MeV on the 4 4 lattice to T tri w 137MeV on the 6 4 lattice. 

Also shown in Figure ||(a) are points from simulations with quark mass am = 0.03. The phase 
boundary turned out to be — within our statistical uncertainties — independent of the the mass. 
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5. The quark number density 

To reveal the properties of the new phase located in the lower right corner of the phase dia- 
gram, we calculated the quark number density, at constant coupling j8 and at constant temperature 
respectively. To obtain the density n q we perform the following integration 



dlndetM \ f , /dlndetM 

dx \-w^ cos{6)) p{x) 



d(«M) / 



The thermodynamic quantity n q are given as usual by 

1 / dlndetM \ 



a 3 N?N t \ d(a/i) / 



(5.2) 



In Figure |2](b) we show the baryon number density, which is related to the quark number density 
by ris = n q /3. The results are plotted in physical units and correspond to a constant temperature of 
T m 143 MeV (4 4 lattice), T » 124 MeV (6 4 lattice) and T « 93 MeV (6 4 x 8 lattice). In order to 
divide out the leading order cut-off effect, we multiply we have multiplied the data with the factor 
c = SB(N t ) I SB, which is the Stefan-Boltzmann value of a free lattice gas of quarks at a given 
value of N t , divided by its continuum Stefan-Boltzmann value. At the same value of the chemical 
potential where we find also a peak in the susceptibility of the plaquette (jU c ), we see a sudden rise 
in the baryon number density. Thus for fi > \i c we enter a phase of dense matter. The transition 
occurs at a density of (2 — 3) x n^, where denotes nuclear matter density. Above the transition, 
the density reaches values of (10 — 20) x n^. Quite similar results have been obtained recently by 
simulations in the canonical ensemble [Eh. 
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